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Abstract 

Molecular dynamics simulations are used to model the thermal properties of a fluid containing 
solid nanoparticles (nanofluid). The flexibility of molecular simulation allows us to consider the 
effects of particle mass, particle-particle and particle-fluid interaction and that of the spatial 
distribution of the particles on the thermal conductivity. We show that the heat conductivity of 
a well dispersed nanofluid is well described by the classical Maxwell Garnett equation model. In 
the case of particle clustering and strong inter particle interactions the conductivity can be again 
described by effective medium calculation taking into account the aspect ratio of the cluster. 
Heat transfer is increased when particles are aligned in the direction of the temperature gradient. 
This kind of collective effects could be a first step to understand the substantial increase in the 
conductivity observed in some experiments. 
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Many experimental studies have suggested that the thermal conductivity of colloidal 
suspensions referred to as "nanoffuids" is unusually high^. Predictions of effective 
medium theories are accurate in some cases^ but generally fail to account for the large 
enhancement in conductivity. In spite of a large number of - sometimes conflicting or 
controversial - suggestions and experimental findings 4 *, the microscopic mechanisms for 
such an increase remain unclear. One of the possibilities that was suggested was the 
effect of Brownian motion 5 , and appeared to be an attractive and generic explanation. 
The essential idea is that the Brownian velocity of the suspended particle induces a fluc- 
tuating hydrodynamic flow^, which on average influences (increases) thermal transport. 
This mechanism is different from transport of heat through center of mass diffusion, the 
contribution of which was previously shown to be negligible^. However, some recent exper- 
imental high precision studies reported a normal conductivity in nanoparticle suspensions 
at very small volume fractions below 1%^, questioning the validity of this assumption. 
Recent simulations also showed that normal conductivity is expected for low volume 
fractions (around 3.3%) 11 as well as for volume fractions up to 13% 10 establishing that 
the physical parameter determining thermal properties should be the particle interfacial 
thermal resistance. 

We also recently established by simulations that Brownian motion and Brownian ve- 
locity field have no effect on the heat transfer of a single nanoscopic particle with the 
surrounding fluid 1 ^. Thus the explanation of the increase in conductivity is currently to 
be looked for mostly in some collective effects between particles - a field that has not been 
studied through microscopic simulation. 

In this work we use non equilibrium molecular dynamics "experiments" to explore 
further the transfer of heat in a model fluid containing nanoparticles. We make use 
of the flexibility allowed by molecular simulations to explore extreme cases in terms e.g. 
particle/fluid mass density mismatch. We concentrate on model systems that are expected 
to be representative of generic properties. We explore a large range of parameters and 
make a quantitative comparison with effective medium calculations. By studying the 
thermal conductivity in a system with two particles and by precisely controlling their 
positions we are able to observe the influence of collective effects consisting in different 
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particle-particle interactions and displacement with respect to the temperature gradient. 

We start by describing our simulation methodology and presenting results about the 
conductivity of a fluid containing a single nanoparticle in a temperature gradient. We 
then study the conductivity of a system containing two particles varying their positions 
and the intensity of their interactions. 



I. SIMULATION METHOD 



The model fluid used in this study is a simple Lennard- Jones liquid. The nanoparticles 
(solid phase) are obtained by a spherical cut of a bulk FCC crystal. All atoms in our 
system interact through Lennard- Jones interactions 

(4e((a/r) 12 -c(c/r)% r < r c 
Uij(r) = { (1) 



1 



0, r > r c 

where r c = 2.5a. The coefficient c is equal to 1 for atoms belonging to the same phase, but 
can be adjusted to modify the wetting properties of the liquid on the solid particle^ 1 ^. 
A coefficient of c = 1 defines a wetting interaction while the non wetting case is modelled 
by c = 0.5. Within the solid particles, atoms are linked with their neighbors through a 
FENE (Finite extension non-linear elastic) bonding potential: 

U FENE (r) = ^RoMl-(^-) 2 ), r<R (2) 

l Ho 

where Rq = 1.5a and k = 30.0e/a 2 . This potential, combined with the Lennard- Jones 
interaction results in a narrow distribution of the distance between linked atoms around 
0.97cr. A solid particle in the fluid was prepared as follows: starting from a FCC bulk 
arrangement of atoms at zero temperature, the atoms within a sphere were linked to their 
first neighbors by the FENE bond. Then the system was equilibrated in a constant NVE 
ensemble with energy value corresponding to a temperature T — 1. A particle contains 
555 atoms, surrounded by the atoms of the liquid. The number density in the system is 
p = 0.85o~~ 3 . As the simulated particles are not exactly spherical, but present some FCC 
facets, their radius was estimated from the radius of gyration: 

<^> = ^X>"'- C «) 2 = ^ (3) 



where R 2 is the measured radius of gyration of the particle atoms, and the second equality 
applies to an ideal sphere. Taking a = 0.3nm this corresponds to a particle radius of order 
Rpart ~ l.bnm. The solid particles obtained by this procedure behave very closely as 
ideal harmonic solids. Non linear effects associated with the non linearity of the bonding 
potential are absent in a wide temperature range. This assumption was verified through 
equilibrium simulations at different temperatures, monitoring the energy per particle. The 
observed relation is linear and indicates a particle heat capacity very close to 3ksTN, 
as for an harmonic ideal solid, in a temperature range from T = 1 to T = 3.5. We 
study systems of wetting ("hydrophilic", c = 1) and non wetting ("hydrophobic" c = 0.5) 
particles, different values of the density mismatch between the solid and liquid phase 
(m p /mi = 1, 50, 100 where mi is the mass of a fluid atom and m p is the mass of a particle 
atom) and different particle volume fractions. 

A. Simulating Heat Flow and Measuring Conductivity 

The most simple and direct method of measuring thermal conductivities in a simulation 
is undoubtedly non equilibrium molecular dynamics. With this set up a temperature 
gradient is established through the sample and thermal properties are calculated from the 
measurement of energy fluxes. This method involves locally injecting and evacuating heat 
in the system. In our study this is achieved by applying two thermostats with different 
temperatures in the two ends of a fluid slab. The periodicity of the box is maintained in 
the directions perpendicular to the temperature gradient and the system is non periodic 
in the direction of the heat flux. If one allows the gradient to change sign in the simulation 
box a stationary heat flow can be achieved in a fully periodic box 11 . However with this 
setup special attention should be payed in the regions of abrupt change in the temperature 
gradient, and the simulated system is less realistic. 

The systems we study are periodic in the x and y directions. In the transverse z 
direction, the liquid is confined by a repulsive potential (ideal flat wall). The thermostats 
are applied to a fluid slice in the vicinity of the walls of width of around 3 atom diameters. 
The thermostat consists in a rescaling of the velocities of the particles currently present 
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in the slice at a given time interval. The temperature measured locally in slices parallel 
to the heat flow direction shows that a linear profile is established in the fluid slab. The 
value of the temperature gradient depends on the temperature of the thermostats and 
their rescaling time constant (see fig. [1]). 




FIG. 1: Dependence of the temperature profile on the thermostat time constant for a pure 
liquid. A smaller rescaling interval increases the efficiency of the energy transfer between the 
thermostat and the liquid inducing a smaller temperature jump. 

A smaller rescaling interval increases slightly the energy transfer between the thermo- 
stat and the fluid thus maintaining the thermostatted layer at an average temperature 
closer to the temperature of the thermostat. It is thus important to have the same time 
constant for the thermostat in all the systems in order to compare the values of the energy 
flow and detect its variations in the different setups. We set this value to dt — 1.5r in the 
simulations. Using this heat transfer setup, the local thermal conductivity is given by: 

a « - M) (4) 

where J(z) is the heat flow and T(z) the temperature. This formula is easily generalized 
to the whole slab by taking the heat flow at the thermostats and a mean slope of the 
temperature profile. The problem is that such an estimation of the mean gradient is 
prone to error. The mean temperature profile for a pure liquid has a well defined slope, 
even if the fluctuations and the fitting procedure cause error. The situation is worse in 
the presence of nanoparticles. The temperature profile averaged over liquid and solid 
atoms has noticeable fluctuations around the particles as the gradient is different in the 



two phases. That is why measuring an effective value of the conductivity that does not 
involve fitting and assuming linear temperature profile is preferable. If the thermostat 
relaxation time is constant for all systems, a well defined values are the two thermostats 
temperatures. Thus we define the slab conductivity as: 

Kff = m - t 2 )/l (5) 

where J is the mean stationary heat flux measured by the thermostats, Tj is the temper- 
ature of the thermostat i and L is the length of the simulation box in the z dimension. 
The conductivity defined in this way is sensitive only to variations in the energy flux and 
is a precise method for capturing variations in conductivities between systems. 

To avoid any effect of thermophoresis or coupling of the thermostat to the particle, 
the particles are constrained to stay away from the thermostatted regions by tethering 
weakly their center of mass to a fixed point by harmonic springs of stiffness k = 30. 
Controlling the particles position also allows us to study different configurations and ex- 
plore the influence of the spatial distribution of the particles on the thermal properties. 
The temperatures of the two thermostats were T\ = 2 and T 2 = 1. In order to compare 
the conductivity results for the different systems they were first equilibrated to the same 
pressure at a temperature T = 1.5, then a non equilibrium run was performed for about 
1500-2000 tlj to make sure the pressure stays the same for the systems of different nature 
and finally a production run of about 15000 t^j during which thermostats energy, particle 
positions and temperature profiles are monitored. In nanoscale systems it has been ob- 
served that interfacial effects are very important^. In order to compare quantitatively the 
conductivity variations to the predictions of effective medium calculations it is necessary 
to know the value of the interfacial (Kapitza) thermal resistance. The values used in our 
study were determined for different wetting and particle masses by a transient adsorption 
simulation as explained in ref.— . In real units, a value Rk = 1 corresponds typically to 
an interfacial conductance G = 1/Rr, of the order of 100MW/Km 2 1 . We study systems 

1 The conversion to physical units is made by taking a Lennard- Jones time unit tlj = 10 _12 s, and 
a length unit a = 0.3nm. The unit for G is energy /temperature/ (length) 2 /time. As the en- 
ergy/temperature ratio is given by the Boltzmann constant fcs, we end up with a unit for G equal to 
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with particle volume fraction of 13% or 12%. The volume fraction is defined as the volume 
of the particle divided by the volume of the fluid outside the thermostats. 

II. RESULTS AND DISCUSSION 
A. Single Particle in a Heat Flux 

First we investigate the effect of the presence of a single nanoparticle on the thermal 
conductivity of the fluid. In this system the center of mass of the particle is held at equal 
distance from the two thermostats (fig. W)- 




FIG. 2: Snapshot of the system used to evaluate the thermal conductivity with a particle of 
13% volume fraction. 

We performed simulations for a small volume fraction ($ ~ 2%), where we were not 
able to detect any change in thermal conductivity compared to the bulk fluid. At a 
higher volume fraction (<3> ~ 13%), on the other hand, we observe a clear decrease in the 
heat conductivity associated with the presence of the nanoparticle (fig. [3]). The Kapitza 
resistance Rk for the considered particles ranges from 1 to 7 so that the associated 

k B /<J 2 /T LJ ~= lO R W/m 2 /s 
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characteristic Kapitza length is in all cases of the order of the particles diameter. This 
means that the decrease must be interpreted in terms of interfacial effects. To quantify 
these effects, we use the Maxwell- Garnett approximation for spherical particles, modified 
to account account for the Kapitza resistance at the boundary between the two media. 
The resulting expression for the effective conductivity^ (see appendix) is 



where A; and X p are the liquid and particle conductivities, $ is the particle volume fraction 
and a = is the ratio between the Kapitza length (equivalent thermal thickness of 
the interface) and the particle radius. This model predicts an increase in the effective 
conductivity for a > 1 and a decrease for a < 1, regardless of the value of the conductivity 
of the particles or of the volume fraction. The prediction depends very weakly on the ratio 
Xp/Xi, less than 1% for 10 < X p /Xi < 100. The minimum value of A e /// A/, obtained when 
a — > oo, is w hile the maximum possible enhancement (for A p — * oo and Rk — ^ 0) is 

YZ§ ■ The Kapitza resistance can be modified by tuning either the liquid solid interaction 
coefficient c, or the mass density of the solid, or a combination of these two parameters. 
Figure [3] illustrates the variation of the measured effective conductivity for several values of 
the Kapitza resistance, determined independently for various values of these parameters. 
It is seen that the observed variation (decrease in our case) in the effective conductivity 
is very well described by the Maxwell-Garnett expression. This expression also allows us 
to understand why the heat conductivity does not vary in a perceptible manner for small 
volume fractions (~ 2%), for which the predicted change would be less than 2%, within 
our statistical accuracy. 

B. Study of Collective Effects 

Next we study a system of volume fraction 12% containing two nanoparticles in various 
configurations, in order to investigate the influence of collective effects. In order to model 
the influence of microscopic particle clustering we study two particles tethered by soft 
springs and forced to stay in "contact", directly interacting with each other (fig. HI We 




(£(l + 2a) + 2)+2*(J(l-a)-l) 



(6) 



(£(1 + 2a) + 2) -*(£(!-«) -1) 
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FIG. 3: Comparison between the ratio of the effective conductivity to the conductivity of the 
pure liquid of the simulated systems and the values obtained from the Maxwell Garnett equation. 

modify the position of the centers of mass of the particles with respect to the temperature 
gradient - the particles are either aligned parallel or perpendicular to the gradient. In 
these two situations the solid phase has an aspect ratio of either a = 2 (parallel) or 
a = 1/2 (perpendicular). We also varied the Lennard- Jones interaction intensity between 
the two particles (e pp ) in order to increase the rate of particle-particle energy transfer. The 
particles-liquid thermal resistance is also modified by choosing different wetting properties 
or different ratios of the masses of the atoms belonging to the two phases. Finally, to 
explore a broader range of configurations, we also modified the distance between the 
centers of mass of the particles. In most configurations the particles are directly interacting 
with each other and in two setups they are separated along the z direction by a layer of 
atoms from the liquid phase, with a thickness of the order of a particle diameter . An 
overview of the systems investigated is shown in table III BL 

First we study systems with aspect ratio a = 2 and strong particle-particle interaction, 
as a function of the interfacial resistance with th efluid(fig. H]). The strong interaction 
is essentially equivalent to a chemical bonding between the particles. Te two particles 
move as a single rigid body and the heat transfer between them is considerably enhanced. 
We can therefore reasonably compare the obtained conductivities with the values an 
effective medium calculation would predict for an ellipsoid with the same aspect ratio. 
A calculation including interfacial thermal resistance for ellipsoidal particles of the same 
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System name 


Aspect ratio a 


£ vv 


nip /mi 


c 


Rk 


Particle-particle Ar 


a2elmlcl 


2.0 


1.0 


1 


1 


0.8 


2Rp 


a2eWmlcl 


2.0 


10.0 


1 


1 


0.8 


2Rp 


a2elm50cl 


2.0 


1.0 


50 


1 


4.0 


2Rp 


a2el0m50cl 


2.0 


10.0 


50 


1 


4.0 


2R P 


a2el0mlc0.5 


2.0 


1.0 


1 


0.5 


3.2 


2R P 


aO.belmlcl 


0.5 


1.0 


1 


1 


0.8 


2R P 


a0.5el0mlcl 


0.5 


10.0 


1 


1 


0.8 


2Rp 


o0.5elm50cl 


0.5 


1.0 


50 


1 


4.0 


2Rp 


a0.5el0m50cl 


0.5 


10.0 


50 


1 


4.0 


2Rp 


elmlcl 


X 


1.0 


1 


1 


0.8 


4Rp 


elm50cl 


X 


1.0 


50 


1 


4.0 


4R P 



TABLE I: Description of the systems containing two particles . The last column indicates the 
distance between the centers of mass of the two particles. 

size, aligned with the thermal gradient, relevant for this situation can be derived from a 
general expression in ref.— . The ratio of the effective conductivity to the conductivity of 
the pure liquid is given by: 



\i 1 - §f3L z 

where $ is the volume fraction, L zz is given by 

Lyy = 1 — 



2{a 2 - 1) 2(a 2 - If I 2 
with aspect ratio a > 1. The parameter (3 is given by 



cosh 1 a 



where 



h + L ZZ (X C - A;) 
X p 

l + £(2 + l/a)^L„ 



(7) 



(8) 



(9) 



(10) 



In the above X p and A; are the conductivities of the particles and the liquid, R p is the 
particles radius and Rk - the Kapitza resistance of the particle-liquid interface. The 
relations hold for an aspect ratio a > 1. 



10 



FIG. 4: Snapshot of the system containing two particles aligned with the temperature gradient 
having strong particle-particle interaction. 

As can be seen in figure [5] the measured conductivities are in reasonable agreement 
with the calculation. For the smallest interfacial resistance there is an enhancement of 
the composite system conductivity by 5 to 10%. For larger values of Rk the small radius 
of the particles results in a decrease in the conductivity, due to interfacial effects. 
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Calculation 
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□ m=50, c=1 
o m=1 , c=0.5 
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FIG. 5: Comparison of the measured variation in conductivity in the e = 10, a = 2 systems to 
the prediction of the effective medium calculation with aspect ratio 2 and volume fraction 0.12 
(equation [7]) . 

Next, we study the thermal behavior of the systems with particles aligned with the 
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thermal gradient but much weaker interactions. The particle-particle interaction is taken 
to be neutral, e pp = 1 (configurations a2elm50cl, a2elmlcl in table III Bl . For this weaker 
interaction, we find that the conductivity is typically 4% below the one obtained with the 
strong attractive interaction. This difference can be understood from the temperature 
profiles shown in figure [61 The weaker interaction results in a higher resistivity in the 
"neck" region, at the boundary between the two particles. As the strongly attractive in- 
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FIG. 6: Temperature profiles in the particles in the m = 1 systems in stationary heat flow. The 
temperature profile within the liquid is also shown for one of the systems. For the other systems 
the liquid temperature profile is very similar. Stronger interaction decreases the temperature 
difference between the particles by approximately 50% compared to particles separated by the 
same distance with neutral interaction. 

teraction decreases the particle-particle thermal resistance, the particles temperatures are 
closer for e pp = 10. As the particles are much more conductive than the liquid phase their 
temperature varies very little within their dimension. We observed that the temperature 
difference between the particles decreases by a factor of two when the interaction intensity 
is increased by a factor of ten: AT = 0.22 for e pp = 1 and AT = 0.12 for e pp = 10 in 
the systems of m = 1. In the m = 50 case AT(e pp = 1) = 0.44 and AT(e pp = 10) = 0.2. 
Because of the layers of fluid (of thickness around 5a) remaining in both cases between 
the particles and the thermostats the global conductivity of the slab is still dominated by 
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interfacial effects and its increase is small. 

When the particle-particle distance is increased, so that the particles do not interact 
directly with each other, the conductivity of the sample stays nearly identical (difference 
~ 1% in favor for the system where the particles are closer) to the case where they are 
in contact and with neutral interaction. Knowing that the contact area is very small 
when the particles are close and also given that the interaction is not strong, being equal 
to the interaction with the liquid atoms e pp = e = 1, the inter particle distance in this 
case does not play an important role in the value of the conductivity. As can be seen 
in fig. El whenever e pp = e = 1 the average temperature of the particle is close to the 
average temperature of the liquid at the z coordinate of its center of mass, the particle 
thermalizes with the surrounding fluid. In contrast, when e pp = 10 the inter particle heat 
flux becomes important and the two particles behave more like a single body. 

We now turn to the conductivity of the systems where the particles are in the plane 
situated in the middle of the box and orthogonal to the temperature gradient. The heat 
fluxes measured for such systems of aspect ratio 1/2, with strong and neutral particle- 
particle interaction, were identical within our statistical accuracy (difference less than 
one percent). The conductivities in this case are slightly lower, with about 4.5% (2% for 
m — 1 and 7% for m = 50) than in the case where the particles are aligned with the 
temperature gradient and has neutral interactions (e pp = 1). The difference is further 
increased to ~ 8% (7% for m = 1 and 10% for m = 50) when the a = 0.5 systems are 
compared with the a = 2 and e pp = 10 systems. 

In summary, we showed that the mutual positions of the particles in suspension in the 
fluid has an influence on the thermal conductivity of the system. If the particles interact 
so that clustering occurs in the suspension, the global conductivity of the nanofluid can 
increase to values higher than the one of a pure fluid, even if a well dispersed suspension at 
the same volume fraction has a conductivity lower than the pure system. The alignment 
with the temperature gradient enhances the conductivity and its effect on a microscopic 
level can be predicted by effective medium calculation taking into account the aspect 
ratio of the particle cluster. According to the effective medium prediction the ratio of the 
slab conductivity to the conductivity of the pure liquid grows essentially linearly with the 
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number of aligned particles (or the aspect ratio) (see fig. [7J) when this number is smaller 
than the ratio of the conductivities of the two phases (A p /A;). Hence even a moderate 
clustering (e.g. in strings of 3 to 4 particles) could be sufficient to interpret increases in 
the thermal conductivity compared to standard effective medium predictions. 




5 10 15 20 25 30 35 40 45 50 

number of aligned particles 



FIG. 7: Prediction of the effective medium calculation (equation [7J for the ratio of the system 
conductivity to the conductivity of the liquid for an interfacial resistance Rk = 1, a volume 
fraction $ = 0.1, and a ratio of 10 between the thermal conductivities of the solid and of the 
liquid. . 

The enhancement tends to (1 — $) + $A p /% when the aspect ratio tends to infinity. 
Thus, following this calculation that seems to give reasonable results at least on the scale 
of two particles and given that the ratio X p / Xi is usually high in real nanofluids, a large 
enhancement can be expected if aggregates of particles are formed in the suspension. Our 
calculation concerns string like aggregates, but it is likely that similar effects would be 
observed with other types of aggregates, e.g. fractal ones with a low enough dimension. 

III. CONCLUSION 

We have explored some aspects of the thermal properties of "nanofluids", at the level 
of model system, individual solid particles and collective effects involving two particles on 
a microscopic scale. 

By varying interaction parameters or mass density, we are able to vary the interfacial 
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resistance between the particle and the fluid in a large range. This allowed us to estimate, 
over a large range of parameters, the effective heat conductivity of a model nanofluid 
in which the particles would be perfectly dispersed. The results for this setup can be 
simply explained in terms of the classical Maxwell- Garnett model, provided the interfacial 
resistance is taken into account. The essential parameter that influences the effective 
conductivity turns out to be the ratio between the Kapitza length and the particle radius, 
and for very small particles a decrease in conductivity compared to bulk fluids is found. 

We also examined the effect of particle clustering and alignment with respect to 
the temperature gradient. We found that alignment improves the conductivity of the 
nanofluid in accordance with calculation with effective medium approach. Increased inter 
particle interaction further enhances conductivity. The results for small clusters of par- 
ticles can again be described by effective medium theory, taking into account anisotropy 
and interface effects. Clustering of particles into string like objects is therefore suggested 
as a possible mechanism for obtaining larger conductivities, compared to the case of well 
dispersed suspensions. 

In order to draw more precise conclusions concerning the conductivity dependence on 
the complex physics of collective effects larger systems containing more complex aggre- 
gates should be examined. The present study provides guidelines and outlines general 
tendencies that are to be expected in a larger and more complex system. 

IV. APPENDIX 

We present here a brief derivation of the Maxwell Garnett equation for the effective 
thermal conductivity in a two phase media (matrix with spherical inclusions) taking into 
account the interfacial thermal resistance. We consider a macroscopically homogeneous 
material with thermal conductivity X e ff in a temperature gradient following some axis: 
T e ff(r) = — g • r. We concentrate on a spherical inclusion of radius r and conductivity 
Ai surrounded by a spherical shell of host material of thickness 7*1, matrix conductivity A2 
and we assume that the inclusion does not change the temperature field for r > r±. The 
two radii define the volume fraction of the inclusion, $ = r^/rf (fig. [8]). The effective 
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FIG. 8: Schematic presentation of the matrix with dispersed inclusions 



conductivity can be determined by solving the steady state diffusion equation for the 
temperature: 

AT(r,#) = (11) 

where 9 is the angle between r and the external gradient g. The solutions of equation [TT] 
in the three different regions are given by: 



Ti(r,0) = Ar cosO, < r < r (12) 
T 2 (r,6) = (Br + E/r 2 )cos6, r < r < n (13) 
T e ff(r,9) = —grcosO, n < r (14) 

The unknown constants in the equations above are to be determined from the appropriate 
boundary conditions. These are obtained by expressing the continuity of the heat flow at 
the domain boundaries and the value of the temperature field: 

Ti(r Q ,e)-T 2 (r ,e) = -\i-^(r )R K (15) 

Axfi(ro) = A 2 §V ) (16) 
or or 

T-AruO) = T.„(r u e) (17) 
A 2 f (n) = Kj^n) (18) 
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where Rk is the interfacial thermal resistance responsible for a temperature jump at 
the matrix-inclusion interface. Substituting the solutions for the temperature fields in 
equations [T6lfT8| we end up with the following relations: 

Arl - Brl - E + \ x R K Ar% = (19) 

AiArjj - \ 2 Brl + 2\ 2 E = (20) 

Br\ + E + gr\ = (21) 

\ 2 Br\ - 2X 2 E + Xeffgrf = (22) 

Using this set of equations, one can obtain the Maxwell Garnett equation for the effective 
conductivity with interfacial thermal resistance: 



K ff _ (&(l + 2a) + 2)+2<fr(%(l-a)-l) 
A 2 (£(i + 2a) + 2)-$(&(l-a)-l) 

where a = RkX 2 /tq. 



(23) 
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